Classical fields method for a relativistic interacting Bose gas 



OO 
O 

o 

(N 
o 

Q 



Oh' 



> 
o 

O 

00 

o 



X 



Emilia Witkowska^, Pawel Zin^''^ and Mariusz Gajda"^'^ 
^ Institute of Physics, Polish Academy of Sciences, al.Lotnikow 32/46, 02-668 Warsaw, Poland 
Institute for Theoretical Physics, Warsaw University, ul. Hoza 69, 00-681 Warsaw, Poland 
Soltan Institute for Nuclear Studies, ul. Hoza 69, 00-681 Warsaw, Poland 
Faculty of Mathematics and Natural Sciences, College of Sciences, Cardinal's Stefan Wyszynski University, 

ul. Dewajtis 5, 01-815 Warsaw, Poland 
(Dated: December 1, 2008) 

We formulate a classical fields method for description of relativistic interacting bosonic particles 
at nonzero temperatures. The method relays on the assumption that at low temperatures the Bose 
field can be described by a c-number function. We discuss a very important role of the cut-off 
momentum which divides the field into dominant classical part and a small quantum correction. 
We illustrate the method by studying a thermodynamics of relativistic Bose field which dynamics 
is governed by the Klein-Gordon equation with Xtp'^ term responsible for interactions. 

PACS numbers: 03.65.Pm, 03.50.-z, 03.75.Nt. 



I. INTRODUCTION 

Achievement of a Bose-Einstein condensation in dilute 
gases of alkali atoms initiated very extensive experimen- 
tal and theoretical studies of these systems. In partic- 
ular, a very important progress have been made in de- 
scription of a Bose-Einstein condensate at nonzero tem- 
peratures. Various rigorous methods were developed to 
study thermal properties of Bose-Einstein condensates 
Unfortunately implementation of these methods is 
very difficult in realistic and physically interesting sit- 
uations. Therefore other approximate approaches were 
considered. These are: the classical fields method based 
on the Gross-Pitaevskii equation (or projected Gross- 
Pitaevskii equation) 0, @|; truncated Wigner method 
stochastic Gross-Pitaevskii equation Q or two fluid 
approach based on coupled Gross-Pitaevskii and kinetic 
equations In our opinion the methods based on the 
classical fields approximation are the most fruitful [9|. 

The classical fields method is based on the observation 
that a granular structure of the bosonic field is totally 
dominated by its wave character at temperatures below 
a condensation point. The bosonic field can be then de- 
scribed by a single complex function. This is the essential 
simplification. It is justified only in a region of param- 
eters were quantum fiuctuations are much smaller than 
thermal ones. Evidently, classical fields method does not 
take into account spontaneous processes or the effect of 
quantum depletion. Advantage of the method is related 
to its great simplicity. It allows for studying not only 
thermal equilibrium properties but also nontrivial dy- 
namics far from thermodynamic equilibrium. It can be 
relatively easy generalized to a system of many interact- 
ing Bose fields. Moreover, the method allows for a direct 
physical interpretation and deeper understanding of un- 
derlying physics. The classical fields method has been ap- 
plied to description of a variety of thermal effects such as: 
a shift of the critical temperature due to interactions [lo| . 
dissipative dynamics of vortexes fragmentation of 
ID condensates crystallization of vortex lattices 



decay of doubly quantized vortex superfluity [15l |. 
dynamics of spinor condensates Kosterlitz-Thouless 
transition in 2D systems [l3|, a phase coherence or 
distillation jl^ of a Bose-Einstein condensate. 

Encouraged by the success of the classical fields 
method we reformulate here this method in order to 
describe interacting relativistic Bose gas. We illustrate 
the method by studying thermodynamics of a relativistic 
Bose field describing system of relativistic massive par- 
ticles and antiparticles whose dynamics is governed by 
the Klein-Gordon equation with Xtp^ term responsible 
for interactions. In a case of relativistic system we have 
to deal with additional degrees of freedom - not only 
the field but also its momentum are dynamical variables. 
Moreover, the total number of particles or antiparticles 
is not conserved. Only their difference - the total charge 
- is preserved during the evolution (20l . [2l| . All these 
differences must be taken into account while formulating 
the classical fields method for relativistic case. 

The system of degenerate relativistic ga s of bosons was 
widely studied both in noninterac ting l22l. l23l. 124 l25l. l26l . 
[27I as well as in interacting cases [28l. I29l.l30r I3ll.l32^l33l|. 
Bose Einstein condensation was predicted for such sys- 
tems. The interacting case is much more interesting as 
some properties of the system can be drastically changed 
[s^ l- In this case the condensation can be view as a 
spontaneous symmetry breaking in the AV'^-theory [33 |. 
The thermodynamic properties of the relativistic Bose- 
Einstein condensate were studied using the path-integral 
method (sil . 32;] and more recently using 211 expan- 
sion p3| in the context of pion and kaon condensation. 
This method was used also for dynamical studies j^. In 
fact a classical approximation, very similar to the one we 
want to introduce here, has been used already in studies 
of relaxation and thermalization in ip'^ theory in 1+1 di- 
mensions [IBl- It was found that the classical approach 
cannot describe this processes accurately as the system 
goes towards a classical not a quantum equilibrium. In 
our opinion this is because the problem of dynamical se- 
lection of the cut-off momentum was not treated with a 
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proper care. The studies of the nonrelativistic systems 
show its extremely important role [s^l- We will discuss 
this issue in detail in Sec. IV. 

In this paper we formulate the classical fields method 
for the relativistic ip'^ theory in 3+1 dimensions. Unlike 
in previously studied 1+1 case [36| the system exhibits 
a Bose-Einstein condensation and particle and antiparti- 
cle spectra are different in the region of small momenta. 
In order to extract a particle aspect of the field from 
the classical wave-like solution we first study the system 
analytically at zero temperature using Bogoliubov ap- 
proximation [111, section III. This enables us to identify 
two types of quasiparticles, their energies and charges. 
This is crucial for further implementation. Then in sec- 
tion IV we formulate the classical fields method for the 
relativistic system and illustrate its application to a uni- 
form system of relativistic, charged Bose field at a ther- 
mal equilibrium. We show how to determine important 
physical parameters, i.e. temperature, chemical poten- 
tial or excitation energies. We stress a very important 
role of an appropriate choice of the cut-off momentum. 
Finally in section V, we show how the physical proper- 
ties of the system vary with temperature. In particular 
we show that the system can be viewed as weakly inter- 
acting gas of quasiparticles even very close to the critical 
temperature. Our particular results agree with the ones 
of [3fl,^,l3a]. 



II. THE MODEL 

We consider a relativistic quantum Bose field '^(r,t) 
in a three dimensional cubic box of volume V = with 
periodic boundary conditions. We assume the standard 
'AV'^' form of the interaction term. Then the Lagrangian 
of the system is L = J dPr L with the Lagrangian density 
of the form: 



C = 



1 
2^ 



A t 



A 



2 ; 



(1) 

where V is the three dimensional gradient, dot denotes 
time derivative and m is a mass {m? > 0) . Our units are 
described in Momentum conjugated to the field is: 



n(r,i) 



(2) 



The field and the momentum satisfy bosonic commuta- 
tion relation: 



i'{r,t),U{r\t) ^2im5^{r~r'), (3) 



while other commutators are equal to zero |20l |. 

Corresponding Hamiltonian of the system is H = 
J cPrH with the Hamiltonian density: 

n ^ —( (2mfUU'' + V^-tV* 



The field evolves in time according to the nonlinear Klein- 
Gordon equation (NLKGE): 



(5) 



The Hamiltonian commutes with the charge operator 
Q — J rf'^r Q, where the charge density Q is defined as: 



2m 



(6) 



Therefore, the charge of the system is preserved in the 
evolution. In this paper we concentrate on the situation 
where the total charge is different than zero, hence with- 
out any loss of generality we assume that the charge is 
positive Q > 0. 



III. BOGOLIUBOV TRANSFORMATION 

We assume that a Bose-Einstein condensate is present 
in the system and the lowest energy mode becomes 
macroscopically occupied. In such a case annihilation 
and creation operators of this mode can be substituted 
by c-numbers. This is the famous Bogoliubov approxi- 
mation [Hi- Within this assumption the field operator 
^{r,t) takes the form: 



#(r,i) 



^00 +X^i(k,t)e»''--^ 



(7) 



where the classical field (/)o represents a condensate, /i is a 
positive energy and A(k, t) are Fourier components (with 
wave- vector k) of the field operator ^ (r, t) . 



A. Excitation energies 

The next step of the BogoHubov method rehes on the 
condition that almost all particles are in A; = mode. 
This makes us to neglect all terms in the Hamiltonian 
that describe the interaction of particles outside the con- 
densate with themselves, keeping only terms responsible 
for interactions of excited particles with the condensate. 
Technically it means that only linear terms in ^(k, t) are 
kept in the nonlinear Klein-Gordon equation. By com- 
paring terms with the same time dependence we get two 
equations. The first relates a value of /i to the condensate 
field: 



= TO2+A|0ol^ 



(8) 



while the second determines dynamics of amplitudes 
i(k,i): 

i(k, t) = 2inA{k, t) ~ Lki(k, t) - A02it(_k^ ^ (9) 

where Lk = + A|0op. The conjugated amplitude 
k, i) satisfies analogical equation. Let us note that 
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the chemical potential is larger than the bare mass so 
the system is Bose condensed |32| consistently with our 
assumption 

These two second order equations have nontrivial so- 
lutions if the corresponding determinant vanishes: 

{ul - - ik)(wk + 2M^k - ik) = iMMY ■ (10) 

The above equation has two solutions (for every mode 
k): 

ojf = + 2^2 ^ A|(/)o|^ - 

- V4kV2 + (2M2 + A|0oP)2, (11) 



and 



u;^ = k^ + 2/i2 + A|0oP + 

+ V4kV + (V + A|0op)2. (12) 

In a long wave-length limit (k — > 0) the frequency wi(k) 
is: 



ivi (k) = kvs , 
where v^ is a sound velocity: 



1/2 



(13) 



(14) 



Evidently these are phonon-like excitations typical for 
massless particles. The linear character of the spectrum 
at small momenta is responsible for a superfluity of the 
system. In the long wave-length limit the second solution 
(IT2I) has the form: 



c^2(k) = k 



2 +1) V2(2/i2 + A|0oP) 



(15) 

This dispersion relation is characteristic for massive par- 
ticles. Let us notice that at zero momentum k — we 
have: wi(0) = and ^2(0) = ^2(2^2 ^ A|(/)op). In a 
short wavelength limit (k ^ 00) both branches of the 
dispersion relation are typical for relativistic noninter- 
acting massive bosons: 



t^i(k) 



0.2 (k) = Jk2 



(16) 
(17) 



where the effective mass is ^^e// = fJ-'^ + MM^- The only 
role of interactions is to modify a bare mass of particles. 

The same spectrum of energies wi(k) and u;2(k) has 
been found inpO| where the path-integral method were 
used and in [33| within 2ni/iV expansion. 



B. Eigenmodes 

The results of the previous subsection allow us to write 
the operator A(k, t) in the form: 

i(k, t) = ake'^^i* + 6ke*"^* + c^^e-"^^' + dy^e"^^* , (18) 



where Ok, &k, Ck, rfk are some operators and wi, W2 are 
positive. Let us notice that terms proportional to oq 
and 60 do not depend on time. Such constant terms are 
already included in 0o- Therefore we have: 



ao = 0, 
So = 0. 



(19) 
(20) 



If A{k, t) satisfies equation ^ then the operators 6k and 
Ck are: 



^k = 7i(k)alk' 
Ck = 72(k)clL^, 



where 



71 (k) 



A</>§ 



wf — 2imjJi — Lk 



72 (k) = 



Ml 



W2 + 2/1C1J2 — 



Thus we have: 



(21) 
(22) 

(23) 
(24) 



i(k,i) = ake-™^*+7ialke™'* + 4e™^*+72d 

(25) 

Commutation relations ^ impose additional constrains 
on operators ak, rfk- To explore this fact we introduce 
scaled operators cik and /3k: 



Ok = /k ak , 
rfk = .gk /3lk ' 



(26) 



where, at this stage, /k and gu are some arbitrary func- 
tions. Consistently with (fT9|) fa is equal to zero: 



/o = 0. 



(27) 



According to ^ we assume bosonic commutation rela- 
tions for the scaled operators: 



(28) 
(29) 



what allows to determine functions /k, 5k- In the mo- 
mentum representation the field operator takes the form: 

i(0,i)e+^^* = 50 [$l{t)+j2{0)Mt)) , (30) 

for k 0, and i(k, t) = ii(k, t) + ij(-k, t) with: 

ii(k,Oe+'''* = /kfak(0+7i(k)«Lk(i)) , (31) 



5k(/3lk(i)+72(k)/3k(t)) , (32) 



for k 7^ 0. In the above equations we used following 
notation ak{t) = ake~^'^i('')* and /3lk(*) = Ply^e+^'^^^^'^K 
Functions /k and 5k are: 



l5k 



l-|7i(k)P' 
2 ^ ^k 

l-|72(k)P' 



(33) 
(34) 
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where: 



^i(k)-c.2(k) 



(35) 



Functions /k and are known up to some arbitrary 
phase factors. In what follows we set the phases to zero. 

Equations l|3ip and ((32l) define the relativistic Bogoli- 
ubov transformation. In the noninteracting limit (A 0) 
operators dj^. and create a bare (noninteracting) par- 
ticle and antiparticle respectively. In the presence of in- 
teractions the system consists of two types of quasiparti- 
cles. The operator dj^ creates quasiparticles with energy 
CL'i(k) and the operator creates quasiparticles with 
energy W2 (k) . The first type of excitations are quasipar- 
ticles which in the small momentum limit are character- 
ized by the phonon-like dispersion relation. Excitations 
of the second type are typical relativistic massive parti- 
cles. These quasiparticles are mixtures of bare particles 
and bare antiparticles. Coefficients 71 and 72 determine 
an amount of this mixing. 



C. Hamiltonian and charge 

In the previous sections we have analyzed the equa- 
tion of motion of the field operator in the Bogoliubov 
approximation. Here we want to approximate the Hamil- 
tonian of the system using arguments of the Bogoliubov 
method. Consistently we keep in the Hamiltonian terms 
proportional to the nonzero value of the field and terms 
proportional to amplitudes up to the second order. 
However, this approximated Hamiltonian does not con- 
serve the field's charge Q. In order to control a mean 
value of the charge we introduce the operator: 



F = H -jlQ, 



(36) 



where /i is a chemical potential and plays a role of the La- 
grange multipHer. The BogoHubov transformation l(30|) - 
(|32l) brings the Hamiltonian to the diagonal form pro- 
vided that the Lagrange multiplier jl in (|36|) is equal to 
the frequency ^ of the zero momentum component of the 
field 0: 



F 



E 



^i(k) 



dfcdk 



(37) 

As we see /i is a chemical potential. In l(37j) we skipped 
the irrelevant constant term — A/(4m)|(/)o|''. The diago- 
nal form of F shows again that we deal with the system 
of noninteracting quasiparticles. A number of the first 
quasiparticles with momentum k is obviously: 



iVi(k) 



dj^dk : 



(38) 



while the number of the second type of quasiparticles has 
the form: 



/^2(k) 



The same approximations bring the total charge to the 
form: 



~ 1 00 1 ' + C + V (gi (k)iVi (k) + 92 (k) 7V2 (k) 



where 



91 (k) 



92 (k) 



and 



i. ^ 



2^ Lk -I- ojI 
2 ?" ' 

2/i Lk + ^2 

0^2 ^ ^1 



(40) 

(41) 
(42) 



(43) 



In obtaining equation l(40|) rapidly oscillating terms have 
been neglected. They will not contribute to a value of the 
charge averaged over time. Quasiparticles of both types 
are carrying a definite charge (in terms of time averaged 
quantities). It turns out that qi{k) is positive while (?2(fc) 
is negative. In the limit of large k the charge of the first 
quasiparticles tends to gi ^ 1 and the charge of the sec- 
ond quasiparticles tends to (72 ^ —1. Assuming thermal 
Bose-Einstein distribution of both kind of quasiparticles 
it is easy to check that the sum in (|40|) converges. The 
quantity C in l(40l) is a quantum depletion. It gives the 
amount of the charge missing from fc = mode in the 
ground state of the many-body system. The quantum 
depletion is proportional to X\(po\'^/m'^ which is a small 
parameter in the case studied here (40l |. Thus at zero 
temperature the whole charge is accumulated in k = 
mode and a zero-momentum field amplitude is: 



(44) 



Equations l(44|) and |[8]) determine values of the chemical 
potential fi and the condensate field intensity |0oP- The 
excitation spectrum and the time evolution of the field 
operator are determined uniquely. 



IV. FORMULATION OF THE CLASSICAL 
FIELDS METHOD 

A. General considerations 



While formulating the classical fields method for rel- 
ativistic bosons we follow the procedure developed for 
nonrelativistic case |;9]. We assume periodic boundary 
conditions and then the natural modes of the system are 
(39) plane waves. Expansion of the field operator and its 
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momentum fl in the basis of plane waves gives: 

^(r,^)=^^'(k,^)e*^ (45) 

k 

fl(r,t) =^^(k,^)e'''^ (46) 

k 

where k is restricted to the first Brillouin zone and takes 
values k = (27r/L)n with rii = 0, ±1, ±2 ... {i — x,y, z), 
and L is a length of the cubic box. If a population of 
a given mode is greater than its quantum fluctuations 
it is legitimate to replace corresponding operators by c- 
numbers. In practice even the mode whose occupation is 
of the order of one can be treated as classical. At zero 
temperature it is only the zero momentum mode k = 
which fulfills this condition. The field operator can be di- 
vided into a large 'classical' part and a small 'quantum' 
corrections. This is the essence of the famous Bogoliubov 
approximation Evidently, at higher temperatures 

many modes are macroscopically populated. Hence, the 
'classical' part of the field includes not only the zero mo- 
mentum component but also components of larger mo- 
menta up to some value k^ax- Occupations of remaining 
modes are set to zero in the lowest order of approxima- 
tion. Consistently, the field operator ^ and correspond- 
ing momentum 11 are approximated by: 

*(r,i)^vl;(r,t)= ^ *(k,t)e^•'^ (47) 

|k|<fc„„, 

n(r,t)^n(r,t)= ^(k,^)e^'^^ (48) 

k|<fe„,„^ 

where ^' and 11 are c-number functions and kmax is a 
cut-off momentum. 

The complex function ^'(r, t) satisfies the same nonlin- 
ear Klein-Gordon equation |[5]) as the field operator. It is 
convenient to rescale the field in such a way that its total 
charge is equal to one, 5" — > ^/Q'^ . After such scaHng 
the equation satisfied by the classical field is: 

\^(r, t) = VH{v, t) - m2*(r, t) - Ag|'J'(r, <)P'J'(r, t) , 

(49) 

and the field ^'(r, i) has the unit charge. Let us no- 
tice that the charge Q and the interaction strength A 
enter the dynamical equation through the product Q\ 
only. However, physical properties of the system at finite 
temperature depend separately on Q and A. These very 
important issue will be discussed later. 

A steady state evolution of the field whose dynamics is 
governed by Eq.(|49|) corresponds to a single realization of 
the relativistic system at a thermal equilibrium. Based 
on the ergodic hypothesis, time averages of physical ob- 
servables with respect to the field fiuctuating around its 
mean value correspond to the microcannonical descrip- 
tion. The only problem is to find a steady state solu- 
tion corresponding to a given energy (temperature) and 
charge. Fortunately, due to the nonhnear term, all so- 
lutions of Eq. l(49|) for a given energy and charge tend to 



the same equilibrium state independently on the partic- 
ular choice of initial conditions. We checked this obser- 
vation by solving numerically the dynamical equations 
for different random choice of the initial field ^'(r, 0) and 
momentum n(r, 0). 

Let us note that because a number of classical modes is 
restricted by the cut-off momentum the equation (|49|) is 
defined on a well specified finite spatial grid Ngrid = ■N'^ 
and the cut-off momentum is kmax — Vii^M /L. There- 
fore, while solving Eq. l(49l) we have to choose a number 
of classical modes first. This number must depend on the 
temperature (or the energy) of the system. However, the 
value of cut-off can be precisely known a posteriori when 
an equilibrium state is generated numerically and ther- 
mal occupations are determined. Therefore, the choice of 
the cut-off momentum is based on a smart guess which 
is verified when the full analysis of solution is done. 

The quantitative criterion for an optimal choice of the 
cut-off momentum is directly related to the fundamental 
assumption of the classical fields method - a macroscopic 
occupation of modes. In order to determine occupation 
of modes we have to be able to assign a unique value of a 
number of particles (which does not have to be an inte- 
ger) to every given amplitude of the k-momentum mode, 
|5'(k)p. In our approach the Hnk between a field ampli- 
tude and a relative occupation of k-momentum mode fol- 
lows directly from the Bogoliubov approximation. There- 
fore, the standard zero temperature Bogoliubov approach 
becomes an inherent part of the classical fields method. 
Various studies of the optimal choice of the cut-off mo- 
mentum in the nonrelativistic theory [1, [s^, |4l| indicate 
that a 'field intensity' (i.e. modulus square) of the high- 
est momentum wave defines the unit which corresponds 
to a nonrelativistic single particle. This allows for de- 
termination of all modes occupations provided that their 
intensities are known. Moreover, as shown in [s^l this 
procedure allows for 'quantization' of the classical field 
and ensures that thermal particle distribution resembles 
the Bose statistics. Such quantization procedure we gen- 
eralize to the relativistic case. 

Let us note that the dynamics can lead to states of mo- 
menta larger than kmax because of the nonlinear term. In 
some implementations of the method a projection oper- 
ator in space coordinates is introduced. This operator is 
applied at every time step to ensure that only projected 
part of the field is considered (l^. In our implemen- 
tation, in order to restrict dynamics to classical modes, 
|k| < kmax^ all momenta are restricted to the first Bril- 
louin zone automatically due to periodic boundary con- 
ditions and usage of the Fast Fourier Transform method 
for time propagation. 

All these general considerations will be elaborated in 
details in the next subsection where we study a thermal 
equilibrium of the relativistic massive charged Bose field. 
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FIG. 1: (Color online) Spectrum of amplitudes. The black 
solid lines refer to the numerical simulations, dashed lines re- 
fer to the peaks' positions given by the Bogoliubov-Popov like 
formulas (see text for explanation). For the zero momentum 
mode the position of the highest peak is IJ,{T) while its am- 
plitude is |(;!)o(r)p. The left peak corresponds to i^^(O) (left 
dashed blue line), and the right peak to ojf (0) (right dashed 
blue line). Notation for the others peaks (fc 7^ 0) is the follow- 
ing: cjf (k) - red left line, cj^(k) - red right line, and <jJ^(k) 
- blue left line, and ii;J(k) - blue right line. Parameters are: 
mass m — 60, \Q = 900, energy E = 75 and the number of 
grid points Af^ = 16^. 



B. Numerical implementation 

Starting from a random initial condition we evolve the 
classical field according to Eq. l(49|) . The first important 
observation is that after some time the field reaches a 
steady state: amplitudes of different momenta modes 
fluctuate around some mean values. Then we look into 
details of the equilibrium state by analyzing the spectrum 
of amplitudes: 



where 



**(x,i). 



(50) 



(51) 



Notice, the spectrum G(k, uj) is proportional to the Green 
function JH]. 

In Fig. [1] we plot spectrum of amplitudes l(50l) for 
several values of momentum k. In order to show 'fine 
structure details' we use the logarithmic scale. Posi- 
tions of characteristic peaks ('singularities') marked by 
the vertical dashed lines correspond to excitation ener- 
gies of quasiparticles. A width of each peak results from 
quasiparticles interactions and is related to a finite life- 
time of a quasiparticle in a given mode. For the zero 
momentum mode (fc = 0) one can observe three distinct 
peaks in the spectrum. The amplitude of the central one 
dominates. The averaged position of the highest peak we 
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k 

FIG. 2: (Color online) Excitation energy of quasiparticles. 
Points refer to peaks' positions calculated numerically accord- 
ing to l|52ll . Lines refer to the Bogoliubov-Popov like formulas 
(see text for explanation), (i) empty circles (upper dashed red 
line) Luf{k), (ii) full circles (lower red line) io^Ck), (Hi) empty 
box (upper dashed blue line) ujti^) ^-^d (iv) full box (lower 
blue line) a;^ (k) . Parameters are the same as in figure [T] 



denote by fJ,{T) (red dashed line), and will refer to it as 
to a chemical potential. The amplitude of this peak can 
be associated with \4>o{T)\'^. Positions of the other two 
peaks are: ^^'^(O) = /i(T) -|-W2(0) (right dashed fine) and 
1^2(0) ~ m(^) ~ '^2(0) (left dashed fine). The spectrum 
of k ^ mode is composed of four peaks. Positions of 
these peaks are symmetric with respect to In gen- 

eral characteristic frequencies visible in the spectrum can 
be divided into two branches: 1) ^^'^(k) — fi{T) ± Li;i(k) 
, and 2) w^(k) = /i(T) ± cjslk). 

Central frequencies of peaks we calculated as follows: 



;(k)^ 



(52) 



where S" is a frequency range defined by the width of the 
peak. In Fig. [2] we show results of numerical simulations 
(marked by circles and squares) and formulas (fTTj) and 
(fT2|) with n{T) and |0o(T)P being determined numer- 
ically from the evolution of the k = mode (marked 
by lines). In the nonrelativistic theory this modifica- 
tion of the zero temperature spectrum is known as the 
BogoHubov-Popov relation (BP) iQ. We notice perfect 
agreement between numerical results and the BP-like for- 
mulas. Let us notice that the first type of excitations 
cji(k) (red lines) is characterized by a linear dispersion 
relation for small |k|. Excitations of the second type, 
UJ2 (k) , (blue lines) have particle like quadratic dispersion 
in the small |k| region. 

If we neglect peak's width then the stationary state of 
the field in the momentum space can be approximated 
by: 

*(k,i)e+''^* = At(k)e-"^^^^^' +A^ik)e"^^^^^' 

+ yl+(k)e-^'^=('')*-h^^(k)e''^^('')*,(53) 
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with amplitudes Af2(M) calculated according to: 

\At,ik)\^ = J2\^ik,C0t,2)\^ 



(54) 



where we sum over all frequencies within a range of peaks 
centered at ujf^2- The form l(53|) is similar to the one given 
by the BogoHubov transformation, Eqs. l(3T|) and i(32l) . 

In order to determine a number of quasiparticles in 
every mode we assume that they are related to ampli- 
tudes of the classical field Eq. (|54|) through the Bogoli- 
ubov transformation. Using formulas l(3T1) and (|32| we 
get the following expression for a number of quasiparti- 
cles of the first type: 



ni(k) 



7Vi(k) _ |^+(k)|^ 



(55) 



0.001 



0.0008 




10 20 30 40 50 60 70 80 90 
k 

FIG. 3: Equipartition relation for the first (full circles) and for 
the second (empty circles) type of quasiparticles. Error bars 
refer to fluctuations of relative energy accumulated in modes 
with k = |k|. The parameters are: Q = 1.25 x 10^ Ncut, 
T = 47.5 iVc«t, m = 60, XQ = 900, E = 75 and ^V^ = 16^ 



and the second type: 



?^2(k) 



jV,(k) _ |^-(k)|^ 

Q 



I5k| 



(56) 



property of a thermal equilibrium state of a classical sys- 
tem. The relative energy per mode is related to the scaled 
temperature T (because we scaled the field): 



According to the definitions Eqs. (|33ll34|) the coefficients 
/k and gk depend on the chemical potential and the con- 
densate amplitude at zero temperature. However, at 
nonzero temperatures we shall use the finite tempera- 
tures values of /i(T) and |(/)o(r)p obtained numerically 
from the spectrum of fc = mode. This is a very impor- 
tant assumption of our approach which determines the 
link between number of particles and the classical field 
amplitudes. When the temperature is very low the as- 
sumption is quite obvious as it is based on the Bogoliubov 
approximation. For higher temperatures this assumption 
is a postulate borrowed from the nonrelativistic imple- 
mentation of the classical fields method. As we will show 
below, this way of obtaining a number of quasiparticles 
leads to equipartition of energy which is typical for a 
thermal equilibrium of a classical system. 

Therefore, the relative energy (normalized to the to- 
tal charge) accumulated in a given momentum mode is 
obviously: 



(ei) ~ (£2) ^ T - T/Q, 



(59) 



£i(k) = cji(k)ni(k) , 
for the first kind of quasiparticles, and: 
£2(k) = W2(k)n2(k) , 



(57) 



(58) 



for the second kind. These relative energies £i(k) and 
£2 (k) are shown in figure [3] as a function of momentum 
k. The energies were averaged over all directions of the 
wave-vector k and error bars represent fiuctuations of 
the relative energy distribution for various directions of 
momentum. The main observation is that the averaged 
energy per mode is constant. Moreover, it is the same for 
both type of quasiparticles £i(k) ~ £2(k). This proves 
the equipartition of energy in the system - the major 



where (£1,2) corresponds to the energy per mode averaged 
over all momenta. The scaled temperature is T = 3.8 x 
10"'' ji^ for parameters of the simulation presented in 
figure [3l 

Having determined the scaled temperature we can find 
its absolute value T if the total charge of the system is 
known. To determine the charge we use the condition of 
macroscopic occupation of all classical modes. According 
to the previous study this can be done by requiring that 
the occupation of the highest momentum mode: 



Nr. 



(60) 



is of the order of one Ncut — 1- This conditions assures 
that the occupation of every other mode considered in 
the evolution {k < kmax) is evidently macroscopic. Eval- 
uated numerically averaged relative occupations of the 
highest momentum mode for both type of quasiparticles 
are almost the same niiVLmax) — n2{}i.max)- For the case 
studied in this section 'ni^2i)<-max) = 8 x 10~^. Thus the 
field's charge is equal io Q — Ncutlni^2i}^max) — 1-25 x 
10^ Ncut what gives the temperature T — QT — 47.5 Ncut 
and the nonlinearity coefficient equal to A = 900 /Q — 
0.0072/Ncut. 

Notice that in the classical fields method the absolute 
temperature T and the total charge Q depend on the 
value of the population at the cut-off momentum. Ar- 
guments leading to the formulation of the classical fields 
method limit a value of Ncut to the number of the order 
of one. In the nonrelativistic case cornparison with the 
ideal gas suggests that Ncut ~ 0.7 (37l.l4l|. 
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C. A choice of the cut-ofF parameter 



The last issue we want to address here is the choice 
of the value of the population at the cut-off momentum 
Ncut- In the nonrelativistic case there is no rigorous way 
for choosing the cut-off parameter. Therefore, some as- 
sumptions for estimation of the value of Ncut are in usage. 
Thus, the value of Ncut is determined by assuming that 
the critical temperature calculated within the classical 
fields method is equal to those of the ideal gas [4l|. This 
assumption seems reasonable in the relativistic case. 

The critical temperature can be determined from the 
analysis of the charge Qo assembled in the zero momen- 
tum mode. This quantity is shown in Fig. [4l Points 
refer to the classical fields results while lines represent 
the best fit. A significant fraction of the charge accumu- 
lates in the zero momentum mode below T ~ 160iVc„t 
and tends to one at T = 0. Let us notice that in a 
finite system, as considered here, there is no phase tran- 
sition. The charge accumulated in the zero momentum 
component of the field does not vanish at any finite tem- 
perature. It rapidly decreases and at some characteristic 
temperature a curvature of the function Qo(T) changes 
sign and the charge slowly goes to zero when the tem- 
perature tends to infinity. This feature can be used to 
define a characteristic temperature. In the nonrelativistic 
case (dilute gas) this characteristic temperature becomes 
the critical temperature in the thermodynamic limit [i^ . 
Unfortunately the classical fields methods cannot be ap- 
plied above transition temperature (no macroscopically 
occupied modes exist) and we cannot explore the region 
of slow decay of the charge accumulated in a zero momen- 
tum field. Therefore we can estimate the characteristic 
temperature by fitting a function to the numerical points 
what gives: Tc ~ ITdNcut- The critical temperature of 
the ideal gas is given by — (SQ/m)^/^, [2a|- Thus by 
comparing values of Tc and we obtain Ncut — 0.2 for 
the studied parameters. 

Determination of the characteristic transition temper- 
ature is not unique for the finite system as considered 
here. Another possibility is based on the analysis of the 
fiuctuations 5Qq = {{Qq) — (Qo)^)^^^ of the charge in 
the zero momentum mode. In Fig. \E\ we show fiuctua- 
tions of the charge SQq versus temperature T. The fiuc- 
tuations increase monotonically with the temperature, 
reach a maximum, and then rapidly fall to a very small 
values. The characteristic transition temperature can be 
defined as the temperature corresponding to the max- 
imum of fiuctuations of the charge Qo- According to 
this definition the characteristic transition temperature 
is To ~ 150 Ncut- The temperature To is smaller than 
Tc (thus the value of Ncut cannot be treated as the uni- 
versal number). This difference results from finite size 
of the system. We expect that in the thermodynamic 
limit both these temperatures become equal as it is in 
the nonrelativistic case |46]. 
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FIG. 4: Charge accumulated in the zero momentum mode 
Qo/Q- Points refer to the numerical simulations while line 
represent the best fitted function. Parameters of the simu- 
lations are: m = 60, Q\ = 900, Q = 1.25 x 10^ Ncut, total 
energy and number of grid points are respectively (from the 
lowest temperature): E = 65, A/"^ = 10^; S = 67, 7V^ = 12^; 
E = 75, Af^ = 16^ E = 96, Af^ = 20^; E = 200, Af^ = 30^; 
E = 360, A"^ = 36^ E = 496, A^^ = 40^ 




200 



FIG. 5: Fluctuations of the charge in the zero momentum 
mode 6Q0/Q as a function of the temperature. Points refer 
to the numerical simulations while the line is plotted just 
to guide an eye and estimate a position of the maximum. 
Parameters of the simulations are the same as in figure |4l 



QUASIPARTICLES AT NONZERO 
TEMPERATURE 



Interactions change a character of elementary excita- 
tions in the system. The first quasiparticles in the small 
momentum limit are phonons and therefore their nature 
is different from relativistic massive particles. For the 
quasiparticles of the second kind the role of interactions 
is to effectively modify a bare mass. Within the classical 
fields method we can investigate how both the interac- 
tions and the temperature affect the physical proprties of 
the system. To this aim, we will use Bogoliubov-Popov 
like relations: Bogoliubov formulas determined in the sec- 
tion III with temperature dependent chemical potential 
/i(T) and condensate amplitude \(j)o{T)\'^. This two last 
quantities are crucial for the whole analysis. In Fig. [6]we 
show /i(r) and |0o(r)p calculated numerically (points) 
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FIG. 6: (a) Chemical potential /i/m and (6) condensate am- 
plitude |0oP/Q- Points refer to the numerical simulations. 
Solid lines are Eq. l[6T|l for the chemical potential, and Eq. lfMjl 
for the condensate field intensity. Parameters of the simula- 
tions are the same as in figure ID 



and we compare them to the results of other authors 
[3^ . [33I I (Hues). The value of the chemical potential a 
increases with the temperature as it is predicted by [32|: 



A /T 



' 3 \ TO 



(61) 



where the zero temperature chemical potential ^ results 
from the Bogoliubov approximation (mass m in denom- 
inator results from our definition of the charge On 
the other hand condensate amplitude |(/)o(r)p decreases 
with the temperature in the way given by the formula 



Q 



(62) 



where \(j)o\^ is the BogoHubov approach result. We ob- 
served quantitative agreement between the classical fields 
method and predictions of [s^l and [s^. Having deter- 
mined /i(T) and |0o(2^)P in the whole range of interesting 
temperatures we can look for thermal properties of the 
system. We concentrate on four quantities: the excita- 
tion energy spectra, decay rates, mixing parameters and 
charges of quasiparticles. 

First, in Fig. [7] we show excitation spectra oji{k) and 
uj20^) for two different temperatures. Points refer to clas- 
sical fields results and lines are: Bogoliubov (solid) and 
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FIG. 7: (Color online) Excitation energy a;i(k) and ti;2(k) as 
a function of momentum = |k|. Blue points refer to the 
numerical simulations for Ei = 75, Ti = 'i7.5Ncut, and the 
number of spatial grid points Af^ = 16^; red points refer to 
the numerical results with E2 = 500, T2 = 163.0A'^cut, and the 
number of grids points JV^ — 40^ . Solid line corresponds to 
the Bogoliubov approximation (marked by BG), and dashed 
line to the Bogoliubov-Popov like formulas. 



Bogoliubov-Popov like (dashed) formulas. For the lower 
temperature, Ti, all results are very close to each other. 
The phonon branch in the first type of quasiparticles ex- 
citation spectrum is clearly visible. For the higher tem- 
perature, T2, the phonon branch practically disappears 
because the sound velocity significantly decreases. These 
features are well known in the nonrelativistic limit [3]. 
The excitation spectrum of the second type of quasipar- 
ticles does not have any visible linear part at long wave- 
length. The finite lifetime of excitations are determined 
by decay rates of quasiparticles. They are related to the 
width of the peaks in the excitation spectrum (see Fig. 
[1]) . The finite width indicates a finite lifetime of quasipar- 
ticles and signifies a weak residual interactions between 
quasiparticles. The peak width can be estimated by: 



ri,2(k) = o.5y^(^?,2(k)>-(c^i,2(k))^ 



(63) 



where (•) represents averaging over frequencies within a 
peak area, compare Eq. (|52l) . In Fig. [8] we show decay 
rates of the first, Fi, and the second, type of quasi- 
particles as functions of k. One can see that in the case of 
the first type of quasiparticles decay rates are increasing 
in the small momenta region where the dispersion rela- 
tion is linear. The decay rates weakly depend on momen- 
tum for large k. This weak dependence is characteristic 
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FIG. 8: (Color online) Decay rates as a function of momen- 
tum: ri(fc) for the first type of quasiparticles and r2(fc) for 
the second quasiparticles. Lines are added to guide the eye. 
Numerical simulations results for two different temperatures: 
Tl = 47.5Ncut, El = 75 and AA^ = 16=* and T2 = 16S.0Ncnt, 
E2 = 500 and Af^ = 40^ 
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FIG. 9: (Color online) Coefficients 7i(k) and 72(k) of the Bo- 
goliubov transformation l|23p and Il24p as a function of fc = |k|. 
Points refer to numerical simulations results for two different 
temperatures: Ti = 47.5iVc«t, Ei = 75 and Af^ = 16^ (blue); 
and T2 = leS.ONcnt, E2 = 500 and = 40^ (red). Solid 
line {BG) refers to the Bogoliubov approximation and dashed 
line to the Bogoliubov- Popov like formulas (BP). 



for the second type of quasiparticles in the whole mo- 
menta range. Moreover, the decay rates of the first type 
of quasiparticles are much smaller than the decay rates 
of the quasiparticles of the second type. In addition, the 
decay rates of both quasiparticles types increase with the 
temperature. 

The relativistic BogoHubov transformation diagonaHz- 
ing approximated Hamiltonian 'mixes' particle's creation 
and annihilation operators: annihilation operator of a 
quasiparticle is a superposition of annihilation operator 
of a bare particle with momentum k and creation opera- 
tor of a bare particle with momentum — k. The amount 
of this mixing is given by the parameter 71 (k) for the first 
type of quasiparticles, and 72 (k) for the second type. The 
coefficients of the relativistic Bogoliubov transformation 
can be calculated numerically in the following way: 



71 (k) 

72 (k) 



i^r(k)i 
l^+(k)l 

l-A2+(k)| 
1-42- (k) I 



(64) 
(65) 



where ^-^e amplitudes defined in l[53|) . On the other 
hand they can be calculated from the Bogoliubov-Popov 
formulas (23), (24). In Fig. [9] we show 71 and 72. Points 
refer to results averaged over angles of the momentum k 



for two different temperatures: Ti = Al.bNcut marked by 
blue, and T2 = IQi.QNcut marked by red. The solid fine is 
the result of the Bogoliubov procedure while dashing fine 
corresponds to the Bogoliubov-Popov like relations. The 
results show that in the linear dispersion regime (phonon 
spectrum) the value of coefficient 71 is close to one. The 
excitations have a quasiparticle character. In a region of 
large momenta 71 is small and quasiparticles are equiv- 
alent to bare particles. For the excitations of the second 
kind 72 is very small for all momenta and decreases with 
increasing temperature. This indicates that distinction 
between the second quasiparticle types and bare antipar- 
ticles is practically meaningless for the values of param- 
eters studied here. 

Finally, let us look for effective charge qi (k) and (72 (k) 
of quasiparticles. In Fig. [TO] we show these effective 
charges as functions of momentum k. All quasiparti- 
cles of the second type have practically the same charge 
qi{k) ~ —1 independently on the value of the tempera- 
ture or momentum. The same charge is carried by an an- 
tiparticle in the noninteracting (A = 0) case [l^l- On the 
other hand the charge qi(k) carried by a single quasipar- 
ticle of the first type strongly depends on its momentum. 
It is greater than one qi> 1 for the phonon-Hke quasipar- 
ticles while it is equal to one for the quasiparticles with 
large momentum. Highly exited quasiparticles of the first 
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FIG. 10: (Color online) Effective charges of the first quasi- 
particle types gi(fe) and the second quasiparticle types 52 (fe) 
for temperatures: 1) Ti = 47.5Ncut (blue points) and 2) 
T2 = 163.0Ai'cut (red points). Points refer to equations ifiTjl . 
I|42p with values oi uji, L02, M and |<?!>o|^ determined numeri- 
cally, dashed lines correspond to the Bogoliubov-Popov like 
formulas, and solid line refers to the Bogoliubov result. 



type have the same charge as particles in the noninter- 



acting case. Our simulations show that the quasiparticle 
charges decreases with increasing temperature. 



VI. SUMMARY 

In this paper we presented the classical fields method 
for description of a relativistic interacting Bose gas at 
thermal equilibrium. First, we studied the system at 
zero temperature using Bogoliubov method. We iden- 
tified two quasiparticles types. Then we formulated the 
classical fields method. We postulated that the classical 
field has a structure of the Bogoliubov transformation. 
This allows us to relate fc-momentum amplitude to the 
occupation number of the fc-mode. By showing equiparti- 
tion of energy and the structure of the mixing parameters 
7i_2 we verified positively our postulate. In addition, we 
discussed some technical details important for the cor- 
rect implementation of the method. The classical fields 
method is an approximate one, similarly to the others, 
and allows for qualitative but nonperturbative study of 
thermal properties of the system. In the case of the rel- 
ativistic interacting Bose gas, the classical fields simula- 
tions show that the system forms a weakly interacting 
gas of two quasiparticles types not only at zero but in 
the whole range of temperatures. 
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